!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!>\brief
!! Make available some flow diagnostics.
!!
!! \details
!!
!<----------------------------------------------------------------------
module mod_dgcomp_flowstate

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_mpi_utils, only: &
   mod_mpi_utils_initialized, &
   mpi_allreduce, wp_mpi, mpi_sum

 use mod_base, only: &
   mod_base_initialized, &
   t_base

 use mod_grid, only: &
   mod_grid_initialized, &
   t_grid

 use mod_atm_refstate, only: &
   mod_atm_refstate_initialized, &
   t_atm_refstate_e, atm_ref_e

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_dgcomp_flowstate_constructor, &
   mod_dgcomp_flowstate_destructor,  &
   mod_dgcomp_flowstate_initialized, &
   update_flowstate, &
   mass_flow, kinetic_energy

 private

!-----------------------------------------------------------------------

! Module types and parameters

! Module variables

 ! public members
 !> \f$\displaystyle K =
 !! \frac{1}{2|\Omega|}\int_{\Omega}\frac{U^2}{\rho}\,dx\f$
 real(wp), protected :: kinetic_energy
 !> \f$\displaystyle \underline{M} =
 !! \frac{1}{|\Omega|}\int_{\Omega}\underline{U}\,dx\f$
 real(wp), allocatable, protected :: mass_flow(:)
 logical, protected :: &
   mod_dgcomp_flowstate_initialized = .false.
 ! private members
 integer :: comm !< MPI communicator
 real(wp) :: gvol !< global volume

 character(len=*), parameter :: &
   this_mod_name = 'mod_dgcomp_flowstate'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_dgcomp_flowstate_constructor(grid,icomm)
  type(t_grid), intent(in) :: grid
  integer, intent(in) :: icomm

  integer :: ierr
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized .eqv..false.) .or. &
       (mod_kinds_initialized    .eqv..false.) .or. &
       (mod_mpi_utils_initialized.eqv..false.) .or. &
       (mod_base_initialized     .eqv..false.) .or. &
! We don't test for mod_atm_refstate_initialized to avoid circular
! dependences; atm_ref_e is not needed in the constructor but it is
! used in update_flowstate.
       (mod_grid_initialized     .eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_dgcomp_flowstate_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   allocate(mass_flow(grid%d))
   comm = icomm
   call mpi_allreduce(grid%vol,gvol,1,wp_mpi,mpi_sum,comm,ierr)

   mod_dgcomp_flowstate_initialized = .true.
 end subroutine mod_dgcomp_flowstate_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_dgcomp_flowstate_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_dgcomp_flowstate_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   deallocate(mass_flow)

   mod_dgcomp_flowstate_initialized = .false.
 end subroutine mod_dgcomp_flowstate_destructor

!-----------------------------------------------------------------------
 
 subroutine update_flowstate(uuu,base,grid)
  type(t_grid), intent(in) :: grid
  type(t_base), intent(in) :: base
  real(wp),     intent(in) :: uuu(:,:,:)
 
  integer :: d, ie, ierr
  real(wp) :: wg(base%m), uuug(size(uuu,1),base%m), rho(base%m), &
    e(base%m), u(grid%d,base%m), senb(grid%d+1), recb(grid%d+1)
  character(len=*), parameter :: &
    this_sub_name = 'update_flowstate'

   d = grid%d
   mass_flow = 0.0_wp
   kinetic_energy = 0.0_wp
   elem_do: do ie=1,grid%ne

     ! quad. weights
     wg = grid%e(ie)%det_b * base%wg

     ! evaluate the solution at the quadrature points
     uuug = matmul( uuu(:,:,ie) , base%p )
     rho = uuug(1,:) + atm_ref_e(:,ie)%rho
     e   = uuug(2,:) + atm_ref_e(:,ie)%e
     u   = uuug(2+1:2+grid%d,:)

     mass_flow = mass_flow + matmul(u,wg)

     kinetic_energy = kinetic_energy          &
       + 0.5_wp*dot_product(sum(u**2,1)/rho,wg)

   enddo elem_do

   senb(1:d) = mass_flow
   senb(d+1) = kinetic_energy
   call mpi_allreduce(senb,recb,d+1,wp_mpi,mpi_sum,comm,ierr)
   mass_flow      = recb(1:d)/gvol
   kinetic_energy = recb(d+1)/gvol
 
 end subroutine update_flowstate
 
!-----------------------------------------------------------------------

end module mod_dgcomp_flowstate

